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Abstract  Phase  field  theory  is  developed  for  solids 
undergoing  potentially  large  deformation  and  fracture. 
The  elastic  potential  depends  on  a  finite  measure  of 
elastic  strain.  Surface  energy  associated  with  fracture 
can  be  anisotropic,  enabling  description  of  preferred 
cleavage  planes  in  single  crystals,  or  isotropic,  applica¬ 
ble  to  amorphous  solids  such  as  glass.  Incremental 
solution  of  the  Euler-Lagrange  equations  corresponds 
to  local  minimization  of  an  energy  functional  for  the 
solid,  enabling  prediction  of  equilibrium  crack  mor¬ 
phologies.  Predictions  are  in  close  agreement  with  ana¬ 
lytical  solutions  for  pure  mode  I  or  pure  mode  II  load¬ 
ing,  including  the  driving  force  for  a  crack  to  extend 
from  a  pre-existing  plane  onto  a  misoriented  cleavage 
plane.  In  an  isotropic  matrix,  the  tendency  for  a  crack 
to  penetrate  or  deflect  around  an  inclusion  is  shown 
to  depend  moderately  on  the  ratio  of  elastic  stiffness 
in  matrix  and  inclusion  and  strongly  on  their  ratio  of 
surface  energy.  Cracks  are  attracted  to  (shielded  by) 
inclusions  softer  (stiffer)  than  the  surrounding  matrix. 
The  theory  and  results  apparently  report  the  first  fully 
three-dimensional  implementation  of  phase  field  the¬ 
ory  of  fracture  accounting  for  simultaneous  geometric 
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1  Introduction 

Brittle  solids — which  include  ceramics,  typical  rocks 
and  minerals,  glass,  and  some  metals — demonstrate  a 
tendency  to  fracture  rather  than  deform  plastically  (e.g. , 
by  dislocation  glide)  when  subjected  to  stresses  exceed¬ 
ing  their  elastic  limit.  Brittle  solids  often  demonstrate 
strong  directional  bonding  at  the  atomic  scale  and  a 
relatively  large  ratio  of  shear  to  bulk  modulus  (Gilman 
2003),  though  exceptions  are  possible.  In  crystalline 
solids,  fractures  may  be  transgranular  (i.e.,  cleavage 
Lawn  1968;  Schultz  et  al.  1994)  or  intergranular  (i.e.,  at 
grain  boundaries),  while  in  amorphous  or  glassy  solids, 
fractures  are  independent  of  intrinsic  microstructure 
but  highly  dependent  on  pre-existing  flaws,  especially 
surface  flaws  (Wilshaw  1971).  In  ceramics  and  glass, 
depending  on  loading  conditions,  fractures  may  consist 
of  dominant  crack(s)  or  distributed  micro-cracks  (Lawn 
et  al.  1994).  Performance  of  ceramics  in  structural 
applications  is  most  often  inhibited  by  their  low  frac¬ 
ture  toughness.  Toughness  can  often  be  improved  by 
modifying  the  micro  structure  to  promote  crack  bridg¬ 
ing  (i.e.,  deflection  or  pull-out  at  second-phase  parti¬ 
cles,  grain  boundaries,  or  other  heterogeneities),  thus 
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shielding  stresses  experienced  by  a  propagating  crack 
and  causing  it  to  meander  or  arrest.  Efforts  to  improve 
ductility  of  ceramics  and  ceramic  nanocomposites  in 
this  manner  include  controlled  introduction  of  weak 
interfaces,  coarsening  and  elongating  the  grain  struc¬ 
ture,  and  embedding  the  polycrystal  with  strong  inclu¬ 
sions  (Faber  and  Evans  1983;  Lawn  et  al.  1994;  Ohji  et 
al.  1998). 

Resistance  to  fracture  is  often  the  controlling  factor 
limiting  applications  in  which  a  brittle  material  can  be 
used;  grain  size  and  morphology  can  be  engineered 
to  optimize  effective  ductility,  impact  resistance,  or 
fatigue  life  for  a  given  application  (Lawn  et  al.  1994; 
Wiederhom  1984).  Selection  or  design  of  a  material 
for  such  an  application  requires  understanding  afforded 
by  a  predictive  model  of  fracture.  Computer  simula¬ 
tions  enable  descriptions  of  fracture  in  brittle  solids 
under  complex  loading  conditions  and  for  nonlinear 
and  anisotropic  material  behaviors,  difficult,  if  not 
impossible,  to  address  using  existing  analytical  meth¬ 
ods  (e.g.,  solutions  available  from  linear  elastic  fracture 
mechanics).  Engineering  finite  element  (EE)  simula¬ 
tions  often  invoke  continuum  damage  mechanics  the¬ 
ories,  wherein  the  tangent  stiffness  of  a  material  ele¬ 
ment  degrades  as  “damage”  accumulates.  Conventional 
continuum  damage  mechanics  theories  (Clayton  and 
McDowell  2003,  2004;  Sun  and  Khaleel  2004;  Clay¬ 
ton  2006,  2008)  require  prescription  of  phenomeno¬ 
logical  kinetic  equations  specifying  the  rate  of  damage 
accumulation;  associated  parameters  in  such  equations 
must  be  tuned  to  experiments  similar  to  those  simulated 
numerically.  Furthermore,  solutions  can  depend  on 
mesh  size.  Simple  models  based  on  the  notion  of  theo¬ 
retical  strength  (Gilman  1960;  Clayton  2009,  2010)  can 
provide  insight  into  directionality  of  fracture  resistance 
but  do  not  enable  prediction  of  morphological  crack 
evolution  for  complex  stress  states.  Cohesive  models 
of  fracture  (Xu  and  Needleman  1993;  Espinosa  and 
Zavattieri  2003;  Clayton  2005;  Arias  et  al.  2007;  Foulk 
and  Vogler  2010)  enable  simulation  of  discrete  cracks 
and  branching;  however,  mesh  construction  should  not 
be  arbitrary  since  crack  paths  are  constrained  to  follow 
element  boundaries.  Extended  EE  methods  (Moes  et 
al.  1999)  alleviate  the  latter  problem  by  permitting  dis¬ 
continuities  to  traverse  elements,  but  predictive  physics 
is  compromised  since  the  crack  propagation  direction 
must  be  specified  by  a  user-defined  criterion.  Simi¬ 
larly,  numerical  methods  involving  incremental  crack 
growth  laws  (Kim  et  al.  1996)  require  specification  of 


criteria  for  crack  extension  such  as  maximum  energy 
release  (Nuismer  1975).  Multi-scale  techniques  (Knap 
and  Ortiz  2003;  Zhang  et  al.  2007)  blending  atomic 
and  continuum  theory  appear  promising  for  describ¬ 
ing  problems  wherein  localized  damage  or  defects  are 
contained  within  a  relatively  small  volume  of  the  entire 
body,  but  are  inhibited  by  difficulties  with  linking  dis¬ 
crete  and  continuous  regions  and  attendant  numerical 
complexity;  purely  atomic  methods  (Knap  and  Sier- 
adzki  1999)  are  necessarily  restricted  to  relatively  small 
system  sizes  and  to  short  time  scales  for  dynamic  sim¬ 
ulations. 

Phase  field  theories  of  fracture  (Jin  et  al.  2001 ;  East- 
gate  et  al.  2002;  Del  Piero  et  al.  2007;  Hakim  and 
Karma  2009;  Kuhn  and  Muller  2010;  Abdollahi  and 
Arias  2012;  Alber  2012;  Borden  et  al.  2012;  Hofacker 
and  Miehe  2012;  Spatschek  et  al.  2011;  Voyiadjis  and 
Mozaffari  2013)  enable  prediction  of  complex  fracture 
processes  without  introduction  of  spurious  physics. 
Apart  from  the  usual  elastic  constants,  essential  mate¬ 
rial  parameters  entering  such  theories  can  usually  be 
directly  related  to  surface  energy  and  diffuse  interfacial 
width,  the  latter  associated  with  regularization  inher¬ 
ent  in  the  phase  field  approach  that  renders  solutions 
mesh  size  independent.  In  essence,  fractures  evolve 
naturally  with  applied  loading  as  the  body  seeks  min¬ 
imum  energy  configurations.  Mathematical  analysis 
(Del  Piero  et  al.  2007)  has  demonstrated  aspects  of  F- 
type  convergence  of  certain  variational  models  towards 
sharp  interface  models  of  Griffith  type  as  interfacial 
width  is  reduced. 

This  paper  develops  a  phase  field  theory  for  fracture 
of  nonlinear  elastic  materials  based  on  a  variational 
approach,  with  numerical  implementation  involving 
incremental  minimization  of  a  suitable  free  energy 
functional.  The  present  novel  model,  whose  governing 
(Euler-Lagrange)  equations  are  derived  in  an  analo¬ 
gous  way  to  an  existing  theory  for  deformation  twin¬ 
ning  (Clayton  and  Knap  2011a),  is  similar  to  a  finite 
strain  variational  model  implemented  in  two  dimen¬ 
sions  (Del  Piero  et  al.  2007),  but  is  new  in  its  incorpo¬ 
ration  of  anisotropy  and  fully  three-dimensional  (3D) 
EE  implementation.  Previous  phase  field  approaches 
have  considered  anisotropic  fracture  energy  (Jin  et  al. 
2001;  Hakim  and  Karma  2009),  albeit  in  the  context 
of  linear  elasticity  and  2D  numerical  simulations.  Fur¬ 
thermore,  these  and  most  other  prior  works  (Eastgate  et 
al.  2002;  Kuhn  and  Muller  2010;  Abdollahi  and  Arias 
2012;  Alber  2012;  Borden  et  al.  2012;  Hofacker  and 
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Miehe  2012;  Voyiadjis  and  Mozaffari  2013)  invoked  a 
dynamic  approach,  requiring  specification  of  one  or 
more  parameters  controlling  the  time  scale  of  frac¬ 
ture  kinetics.  In  contrast,  the  present  model  invokes 
incremental  energy  minimization,  enabling  the  solution 
of  quasi- static  fracture  problems  without  the  need  for 
kinetic  parameter(s).  This  work  represents  the  first  fully 
3D  implementation  of  a  phase  field  theory  of  fracture 
accounting  for  geometric  nonlinearity,  nonlinear  elas¬ 
ticity,  and  possible  surface  energy  anisotropy.  Results 
of  numerical  simulations  of  mode  I  and  mode  II  loading 
validate  the  model.  Simulation  results  on  crack  deflec¬ 
tion  around  or  penetration  through  a  spherical  inclu¬ 
sion  then  follow,  providing  new  insight  into  crack- 
inclusion  interactions  in  nonlinear  elastic  solids  not 
available  from  2D  linear  analytical  solutions  (Tamate 
1968;  Atkinson  1972;  Erdogan  et  al.  1974;  Evans 
1974). 


2  Phase  field  theory 

2.1  Governing  equations 

Einite  deformation  theory  (Clayton  2011;  Clayton  and 
Knap  2011a,b,  2013)  is  used,  where  x  =  x(X)  = 
X  -h  u(X)  are  spatial  Cartesian  coordinates  of  a  mate¬ 
rial  particle  initially  at  X.  Let  V  denote  the  material 
gradient  operator.  The  deformation  gradient,  symmet¬ 
ric  finite  deformation  tensor,  and  ratio  of  deformed  to 
initial  volume  are,  respectively, 

F  =  Vx  =  1  +  Vu,  C  =  F^F,  J  =  VdetC.  (1) 

Degrees  of  freedom  are  displacement  u(X)  and  order 
parameter  r](X),  where  r]  =  0  for  perfect  material, 
T]  =  1  for  fully  fractured  material,  and  ^7  G  (0,  1)  within 
diffuse  interfaces  between  intact  and  failed  domains, 
for  example.  The  total  free  energy  functional  for  the 
body  ^2  is 

\]/(u,r])=  [  [W(Vu,  T]) —T]^ 

JQ  I 

+  /r  :  V?7  (g)  V?7]d^2.  (2) 

Elastic  strain  energy  per  unit  reference  volume  is  W 
and  may  degrade  with  increasing  rj,  surface  energy  per 
unit  reference  area  is  F,  constant  I  is  proportional  to 
the  equilibrium  crack  width  or  thickness  of  the  diffuse 
interface,  and  ^  is  a  symmetric  second-order  material 
property  tensor,  in  what  follows  of  the  form 


ic  =  r/[l-h^(l-m(8)m)],  (3) 

with  m  a  unit  normal  vector  to  a  preferred  cleavage 
plane  in  a  crystal,  for  example.  Parameter  p  penalizes 
fractures  on  planes  with  orientations  different  than  m . 
Eor  isotropic  surface  energy,  e.g.,  as  in  a  glass  or  a  rep¬ 
resentation  of  a  macroscopically  homogeneous  brittle 
isotropic  solid,  ^0  =  0.  Erom  functional  (2),  respective 
Euler-Lagrange  equations  in  Q  and  boundary  condi¬ 
tions  on  its  external  surface  are  derived  using  stan¬ 
dard  variational  methods  and  the  divergence  theorem 
(Clayton  and  Knap  2011a): 

V-P  =  0,  dW/dr]  -\-2rr]/l  =  2W  ’KWtj;  (4) 
t  =  Pn,  h  =  2ic  :  Vr]  (g)  n.  (5) 

Eirst  Piola-Kirchhoff  stress  isP  =  dW/dF  = 
dW/dVu  and  is  generally  not  symmetric,  traction  on  a 
surface  element  with  reference  normal  n  is  t,  and  h  is  sl 
conjugate  force  to  r]  on  boundary  9^2  (note  that  h  =  0 
along  a  free  surface). 

2.2  Elasticity 

Similarly  to  previous  phase  field  models  for  twinning 
(Clayton  and  Knap  2011b,  2013),  compressible  neo- 
Hookean  elasticity  is  considered  with  the  following 
strain  energy  density  function: 

W(C,  r))  =  l[;.(ln  jf  -  ix{2  In  /  -  trC  +  3)].  (6) 

Lame  coefficients  (jl  and  A  depend  on  order  parameter 
Y]  and  possibly  volume  change  measure  /: 

(7) 

Shear  modulus  /x  degrades  from  its  reference  value  /xq 
to  a  minimum  value  ^/xq  in  fully  fractured  domains, 
where  0  <  ^  ^  1.  Bulk  modulus  k  degrades  in  tension 
but  not  in  compression,  in  order  to  prohibit  interpene¬ 
tration: 

k  =  ^o{[^  +  (1  -  0(1  -  04{7  -  1>  +  (1  -  7)*)- 

(8) 

The  initial  (not  degraded)  bulk  modulus  is  Aq  =  Aq  + 
|/xo.  The  following  notation  applies:  (x)  =  IVv  > 
0,  (x)  =  OVv  <  0,  (v)*  =  IVv  >  0,  and  (v)*  = 
OVv  <  0.  Poisson’s  ratio  in  the  undamaged  solid  is 
V  =  Ao/ (2Ao +2/xo) .  The  stress  tensor  is  then,  invoking 
(6), 
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P  =  jjiF  +  (A  In  /  —  ijl)F  ^ 

+  ^J(lnjf(dk/dJ)F-'^.  (9) 

The  rightmost  term  accounts  for  tension-compression 
asymmetry  in  the  degradation  of  the  bulk  stiffness. 
Anisotropic  elastic  constants  are  not  included  here  but 
may  be  considered  in  the  future  by  adapting  methods 
in  Clayton  and  Knap  (201  la);  however,  in  many  single 
crystals,  surface  energy  anisotropy  is  thought  to  domi¬ 
nate  elastic  anisotropy  (Lawn  1968). 

For  purposes  of  comparison  with  known  analytical 
solutions,  linear  elasticity  is  also  sometimes  considered 
in  this  work,  whereby  the  strain  energy  density  function 
is 

W(Vu,  11)  =  L(V  ■  h)2  +  l/itr{[(V«  +  (V«)T)]2}, 

(10) 

and  elastic  coefficients  degrade  in  the  same  way  as  in 
(7)  and  (8),  with  volume  change  in  the  linear  approx¬ 
imation  J  ^  1  +  V  •  M.  The  stress  tensor  [derivative 
of  (10)  with  respect  to  Vm]  in  the  linear  theory  is  sym¬ 
metric: 

1  9 

P  =  A(V  •  M)l-h2/x(VM)sym  +  -(V  •  Ur(dk/dWu)l. 

(11) 

Comparison  of  nonlinear  and  linear  results  also  enables 
evaluation  of  possible  importance  of  nonlinear  aspects 
of  the  model  and  thereby  may  suggest  domains  of  valid 
loading  regimes  for  the  linear  theory. 

The  present  model  permits  degradation  of  the  bulk 
modulus  only  when  volume  change  is  tensile  and  degra¬ 
dation  of  the  shear  modulus  regardless  of  whether  load¬ 
ing  is  tensile  or  compressive.  This  approach,  which  is 
analogous  to  that  implemented  in  Amor  et  al.  (2009) 
and  one  described  in  Spatschek  et  al.  (2011),  enables 
the  material  to  lose  shear  strength  even  when  all  three 
principal  strain  components  are  negative.  A  different 
theory  implemented  in  Miehe  et  al.  (2010),  Hof  acker 
and  Miehe  (2012),  Borden  et  al.  (2012)  in  the  small 
strain  setting  decomposes  the  strain  tensor  into  posi¬ 
tive  and  negative  parts  following  diagonalization.  The 
elastic  strain  energy  dependent  on  the  positive  (tensile) 
components  degrades  upon  fracture,  while  the  com¬ 
pressive  energy  is  unaffected  by  damage  (i.e.,  does 
not  vary  with  changing  values  of  the  order  parame¬ 
ter).  In  contrast  to  the  present  approach,  shear  strength 
loss  does  not  occur  when  all  three  principal  strains  are 
compressive  (Borden  et  al.  2012).  Suitability  of  one 


approach  over  the  other  likely  depends  on  the  material 
and  loading  regime,  and  could  be  evaluated  by  con¬ 
sidering  the  failure  behavior  of  a  material  subjected  to 
triaxial  compression  (Murrell  1965)  or  shock  compres¬ 
sion  under  lateral  pre-stress  (Clayton  2014).  In  brittle 
solids  such  as  ceramics  and  rocks,  when  confining  or 
lateral  stress  is  not  too  large,  shear  fractures  may  occur 
under  compression,  but  extreme  lateral  stress/pressure 
may  suppress  propagation  of  mode  II  cracks  due  to  fric¬ 
tional  resistance  that  tends  to  increase  with  increasing 
confinement  (Murrell  1965;  Clayton  2010). 

3  Numerical  implementation 

The  theory  is  implemented  numerically  using  Lagran- 
gian  finite  elements,  following  general  procedures  out¬ 
lined  in  Clayton  and  Knap  (2011a).  Incremental  solu¬ 
tions  to  (4)  are  sought  using  conjugate  gradient  min¬ 
imization  of  total  energy  functional  ^  of  (2),  subject 
to  constraints  associated  with  external  boundary  con¬ 
ditions  on  9^2.  Let  8r](X)  be  a  local  change  in  order 
parameter  r]  at  material  point  X  induced  by  incremental 
loads.  The  additional  internal  constraint  8r](X)  >  0  if 
r](X)  >  0.9  is  used  to  render  fracture  irreversible  (Del 
Piero  et  al.  2007). 

In  the  current  implementation  of  the  phase  field  the¬ 
ory,  cracks  represented  by  positive  values  of  the  order 
parameter  are  predicted  to  follow  paths  dictated  by 
incremental  total  energy  minimization,  subject  to  the 
irreversibility  constraint  described  above.  When  this 
constraint  is  active,  the  incremental  energy  minimiza¬ 
tion  problem  can  be  viewed  as  minimization  of  energy 
of  an  alternative  system  with  time  dependent  boundary 
conditions  associated  with  introduction  of  new  surfaces 
along  which  rj  >  0.9  is  prescribed;  equilibrium  equa¬ 
tions  (4)  remain  satisfied  in  solutions  thus  obtained  for 
this  alternative  system.  If  the  irreversibility  constraint 
on  5?7  is  not  enforced,  then  damage  is  reversible  and 
cracks  will  heal  fully  upon  unloading,  a  feature  noted 
in  other  phase  field  models  (Hakim  and  Karma  2009). 
If  damage  is  regarded  as  completely  irreversible,  then 
dissipated  energy  can  be  defined  as  the  contribution 
to  ^  of  the  final  two  terms  in  the  integrand  of  (2) 
corresponding  to  the  surface  energy  of  fracture  asso¬ 
ciated  with  the  order  parameter  and  its  spatial  gradient. 
This  is  consistent  with  definitions  of  dissipated  energy 
used  elsewhere  (Miehe  et  al.  2010;  Borden  et  al.  2012; 
Spatschek  et  al.  2011)  when  viscous  effects  are  omit- 
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ted.  Another  definition  of  dissipated  energy,  valid  for 
partial  reversibility  of  damage,  could  be  constructed  for 
a  loading-unloading  cycle:  net  dissipated  energy  could 
be  measured  by  the  difference  between  work  expended 
and  then  recovered  by  the  boundary  conditions,  minus 
any  residual  energy  remaining  in  the  system  upon  exter¬ 
nal  unloading.  Further  review  of  capabilities  of  phase 
field  approaches  to  predict  crack  paths  can  be  found  in 
Spatschek  et  al.  (201 1). 

In  order  to  make  the  results  broadly  accessible  and 
applicable  to  a  number  of  classes  of  brittle  materials, 
numerical  analyses  are  performed  later  in  the  context 
of  dimensionless  parameters  rather  than  property  sets 
peculiar  to  any  specific  solid.  Energy  density  can  be 
normalized  by  /xq,  such  that  W/fio  depends  only  on 

V  and  other  terms  in  the  integrand  of  (2)  depend  on 
dimensionless  constants  F  =  F I /jLqI,  and  Tak¬ 
ing  ^  =  0.01  and  I/Rq  =  constant,  where  Rq  is  a 
fixed  characteristic  length  associated  with  the  problem 
of  interest,  dimensionless  solutions  depend  only  the 
choice  of  {y,  r,  where  ^0  =  0  for  isotropy.  Proper¬ 
ties  representative  of  brittle  solids  are  later  chosen  as 

V  =  0.25  and  F  =  0.01,  while  ^  and  m  are  explored 
parametrically. 

In  order  to  evaluate  accuracy  and  validity  of  the 
model,  simulations  on  a  block  of  elastic  material  with 
a  pre-existing  straight  crack/notch  are  performed.  The 
block  is  of  reference  dimensions  50/?o  x  507?o  x  257?o, 
where  Rq  =  21  is  the  finite  radius  of  the  initial  notch 
tip,  and  the  initial  crack/notch  is  of  length  a  =  25  Rq. 
The  plane  of  this  pre-crack  is  X2  =  0  in  local  Cartesian 
reference  coordinates,  and  let  (r,  0)  he  local  reference 
polar  coordinates  with  origin  at  the  notch  tip.  Along  the 
outer  boundary,  displacement  boundary  conditions  are 
imposed  for  either  pure  mode  I  or  pure  mode  II  loading 
(Zhang  et  al.  2007;  Rice  1968);  e.g.,  for  mode  I: 


where  the  load  parameter  is  A  =  K\l2ii^.  Similar 
equations  apply  for  mode  II,  where  A  =  ^ii/2/xV^ 
(Clayton  and  Knap  2013).  For  a  sharp  crack  in  an  infi¬ 
nite  medium,  linear  elastic  fracture  mechanics  predicts 
crack  extension  will  occur  when  A  >  Kc/2/xv^, 
where  Kq  =  \l2\jiGI(\  —  v),  and  G  =  2 C  for  an 
isotropic  material.  Because  of  the  finite  domain  and 


crack  radius,  values  of  Kyu  entering  A  are  corrected 
herein  by  computing  the  ^-integral  (Rice  1968)  numer¬ 
ically.  Linear  elastic  fracture  concepts  (e.g.,  3  =  G) 
are  used  only  to  define  boundary  conditions  and  aid  in 
analysis  and  interpretation  of  results;  the  phase  field 
simulations  do  not  directly  incorporate  linear  elastic 
fracture  mechanics.  To  mimic  plane  strain  1/3  =  0  is 
imposed  along  the  boundary,  and  h  =  Ois  imposed  on 
all  external  surfaces.  Finite  element  meshes  consist  of 
~  4M  tetrahedral  elements  of  size  small  relative  to  I  to 
resolve  order  parameter  gradients  at  fracture  surfaces. 

4  Numerical  results 

4.1  Cracking  for  isotropic  and  anisotropic  surface 
energies 

Propagation  of  a  pure  mode  I  crack  in  a  notched  body 
with  isotropic  surface  energy  =  0)  is  shown  in 
Fig.  la,  b  for  nonlinear  elasticity  for  two  different  load 
increments,  and  in  Fig.  Ic  for  linear  elasticity.  Initia¬ 
tion  at  the  notch  tip  occurs  at  Ki/Kq  ^  1.03  for  both 
nonlinear  and  linear  elasticity,  validating  the  model  for 
isotropic  surface  energy,  and  neo-Hookean  elasticity 
predicts  a  10-15  %  longer  crack  than  linear  elasticity 
at  the  larger  applied  displacement.  Figure  Id  shows  the 
effect  of  anisotropic  surface  energy  with  =  100  for 
a  cleavage  plane  twisted  at  angle  (p  =  jt  14  relative  to 
the  pre-existing  crack  plane;  i.e.,  m  =  [cos  0,  0,  sin  0], 
where  normalization  of  the  loading  for  twist  is  (Wieder- 
horn  1984) 

Kq  =  y2/xG/(l  —  y)  sec^  0  (twist).  (14) 

Figure  le  shows  a  similar  result  for  a  cleavage  plane 
tilted  at  angle  0  =  7r/4,  i.e.,  m  =  [cos0,  sin0,  0], 
where  normalization  for  pure  tilt  misorientation  is  by 
the  factor  (Wiederhom  1984;  Xu  et  al.  2003) 

Kc  =  s/2f^GI{l  -  v)  sec^((t>/2)  (tilt).  (15) 

Shown  in  Fig.  2a,  b  are  components  of  first  Piola- 
Kirchhoff  stress  tensor  P  of  (9),  normalized  by  the  ini¬ 
tial  shear  modulus.  Stress  concentrations  are  of  large 
magnitude  (e.g.,  locally  on  the  order  of  10%  of  the 
modulus)  but  remain  bounded  in  part  due  to  the  finite 
radius  of  the  crack  tip  resolved  numerically. 

As  validated  in  Fig.  3,  the  phase  field  model  cor¬ 
rectly  predicts  that  twist  misorientation  resists  crack 
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(a)  (b) 


Fig.  1  Phase  field  predictions  for  pure  mode  I  or  mode  II  loading, 
with  fractured  material  (r]  >  0.7)  removed  for  visualization:  a 
nonlinear  elastic,  ^  =  0,  Ki/Kq  =  1.1  b  nonlinear  elastic, 
yd  =  0,  Ki/Kc  =  2.4  c  linear  elastic,  yd  =  0,  Ki/Kq  =  2.4 
d  nonlinear  elastic,  0  =  7r/4  (twist),  yd  =  100,  Ki/Kq  =  1.1 
e  nonlinear  elastic,  (j)  =  n  14  (tilt),  yd  =  100,  Ki/Kc  =  1.5  (f) 
linear  elastic,  yd  =  100,  Ku/Kq  =  1.5 


propagation  more  than  the  same  angle  of  tilt  misorien- 
tation  (Wiederhorn  1984).  Figure  If  depicts  a  result  for 
pure  mode  II  loading  with  yd  =  100  set  to  favor  crack 
propagation  on  m  =  [1,  0,  0].  In  all  cases  shown  in 
Figs.  1  and  2,  when  imposed  boundary  displacements 
are  large  enough  that  cracks  become  overdriven  (e.g., 
K  >  they  do  not  necessarily  propagate  com¬ 

pletely  through  the  specimen  in  an  unstable  manner 
because  of  the  boundary  conditions  that  tend  to  repel 
the  crack  as  it  approaches  the  lower  edge  and  possi¬ 
ble  effects  of  nonzero  ^  that  serves  to  maintain  some 
residual  stiffness  in  the  fully  damaged  material. 

In  the  examples  involving  misorientation  (Fig.  Id, 
e),  cracks  demonstrate  non-uniform  fronts  suggestive 


Fig.  2  Phase  field  predictions  for  pure  mode  I  loading,  nonlinear 
elastic,  yd  =  0,  Ki/Kc  =  2.4,  with  fractured  material  (r]  >  0.7) 
removed  for  visualization:  a  normal  stress  component  Pn/l^o  = 
PxX /)^o  b  normal  stress  component  P22/ =  PyY /l^o  [referred 
here  to  global  coordinates  with  (X,  Y)  oriented  (normal,  parallel) 
to  the  crack  face] 


Fig.  3  Predictions  of  phase  field  model  (nonlinear  elastic,  yd  = 
100)  for  extension  of  a  mode  I  crack  on  plane  misoriented  by 
an  angle  0  of  twist  or  tilt  from  the  original  cleavage  surface, 
compared  to  linear  elastic  fracture  mechanics  solution  (Gell  and 
Smith  1967;  Wiederhorn  1984;  Xu  et  al.  2003).  is  the 

ratio  of  driving  force  for  extension  on  the  misoriented  plane  rel¬ 
ative  to  that  for  0  =  0° 

of  possible  contributions  of  mode  III  loading.  Presum¬ 
ably,  such  effects  arise  from  the  finite  thickness  of  the 
domain  in  the  X3  direction,  since  far-field  in-plane 
boundary  conditions  (12)  and  (13)  apply  for  pure  mode 
I  loading  of  the  initial  notched  configuration  under 
plane  strain.  For  the  case  shown  in  Fig  Id,  dependence 
of  the  solution  on  X3  is  expected  regardless  since  the 
cleavage  plane  normal  m  has  a  nonzero  out-of-plane 
component,  though  the  predicted  cusps  in  the  crack 
profile  may  be  affected  by  or  attributed  to  mode  Ill- 
type  contributions.  Analytical  solutions  to  which  the 
simulation  results  are  compared  in  Fig.  3,  as  stated 
in  Wiederhorn  (1984)  and  derived  in  Gell  and  Smith 
(1967),  Xu  et  al.  (2003)  do  not  account  for  finite  thick¬ 
ness  effects.  Their  agreement  with  initiation  predicted 
by  the  phase  field  simulations  is  reasonable,  however, 
and  it  is  noted  that  non-uniformity  appears  less  pro¬ 
nounced  at  the  onset  of  crack  growth.  The  present 
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phase  field  theory,  like  the  cited  analytical  solutions, 
prescribes  a  single  fracture  energy  for  a  given  material 
or  interface,  and  does  not  distinguish  among  possibly 
different  work  of  separation  (toughness)  values  under 
modes  I,  II,  and/or  III.  In  contrast,  some  cohesive  mod¬ 
els  of  fracture  such  as  those  described  in  Espinosa  and 
Zavattieri  (2003)  permit  the  freedom  of  different  val¬ 
ues  of  work  of  separation  under  different  pure  modes 
of  loading  as  well  as  under  mixed  mode  conditions. 


4.2  Crack  interaction  with  a  spherical  inclusion 

A  second  set  of  phase  field  simulations  considers  the 
presence  of  an  elastic  inclusion  embedded  in  an  other¬ 
wise  homogeneous  matrix  material  of  the  same  geome¬ 
try  (i.e.,  with  notch  or  pre-crack  of  length  <2  =  50/),  and 
subjected  to  the  same  mode  I  displacement  boundary 
conditions  as  considered  in  Figs.  1  and  2.  The  spheri¬ 
cal  inclusion  has  initial  radius  Ri,  is  centered  at  a  dis¬ 
tance  25/  from  the  initial  notch  tip,  and  here  obeys  the 
same  general  constitutive  laws  as  the  matrix  material, 
though  inclusion  and  matrix  may  have  different  elas¬ 
tic  stiffness  and  different  dimensionless  fracture  energy 
r.  Both  matrix  and  inclusion  are  assigned  v  =  0.25 
(i.e.,  Xq  =  /xo)  and  =  0  (isotropy).  Depending  on  the 
mismatch  in  properties  of  the  two  materials,  a  propagat¬ 
ing  mode  I  crack  will  either  penetrate  the  inclusion  or 
deflect  around  it.  Representative  examples  of  the  latter 
are  shown  in  Fig.  4  (order  parameter  r])  and  Fig.  5  (nor¬ 
mal  and  shear  stress  components  for  a  planar  slice  of  the 
three-dimensional  domain).  As  demonstrated  in  Fig.  5, 
both  normal  and  shear  stresses  are  large  but  bounded 
near  the  crack  tip  as  it  propagates  around  the  stiff  and 
strong  inclusion.  Recall,  on  the  other  hand,  that  in  lin¬ 
ear  elastic  fracture  mechanics,  singularities  arise  at  the 
crack  tip  (with  shear  stress  vanishing  atO  =0)  for  pure 
mode  I  loading  (Rice  1968).  The  present  phase  field 
simulations,  which  omit  anisotropy  of  fracture  tough¬ 
ness  and  elastic  constants,  are  physically  representa¬ 
tive  of  a  brittle  glassy  material  with  perfect  bonding 
between  inclusion  and  matrix,  for  example. 

Analytical  solutions  for  crack-inclusion  interac¬ 
tions  (Tamate  1968;  Atkinson  1972;  Erdogan  et  al. 
1974)  idealize  the  problem  as  2D  and  linear  elastic,  but 
do  enable  prediction  of  effects  of  elastic  properties  and 
geometry  on  local  stress  intensity  factors.  These  works 
all  predict  a  decrease  in  crack  driving  force  as  the  stiff¬ 
ness  of  the  inclusion  increases.  A  similar  effect  is  pre- 


Fig.4  Phase  field  prediction  of  crack  defiection  around  a  strong 
second-phase  inclusion  of  radius  Ri  =  1.51  for  far-field  mode  I 
loading  (r/Po  =  5,  k/ko  =  1)  with  failed  material  (rj  >  0.7) 
removed  to  visualize  crack  propagation 


dieted  by  the  phase  field  model,  as  shown  in  Fig.  6:  for 
the  same  far-field  boundary  conditions,  the  crack  moves 
closer  to  the  soft  inclusion  than  the  stiff  one.  The  present 
phase  field  approach  predicts  evolving  crack  geometry 
(i.e.,  growth)  as  the  specimen  deforms  to  possibly  large 
local  strains,  whereas  analytical  results  (Tamate  1968; 
Atkinson  1972;  Erdogan  et  al.  1974)  only  indicate  a 
tendency  for  crack  extension  from  a  pre-existing  flaw 
of  fixed  geometry. 

Numerous  simulations  of  this  inclusion-crack  inter¬ 
action  problem  are  performed,  wherein  elastic  stiffness 
k  and  fracture  energy  F  of  the  inclusion  are  varied,  with 
analogous  properties  Xq  =  /xq  and  Fq  of  the  matrix 
held  fixed,  and  with  the  same  far-field  mode  I  bound¬ 
ary  conditions  imposed  in  each  case.  Results  enable 
quantification  of  effects  of  properties  on  propagation 
behavior  (i.e.,  defiection/bridging  versus  particle  frac- 
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Fig.  5  Phase  field  predictions  for  stresses  along  mid-plane  X3  = 
0  for  crack  deflection  around  a  strong  and  stiff  second-phase 
inclusion  {R\/l  =  F / Fq  =  k/ko  =  5)  with  failed  material 
{rj  >  0.7)  removed  to  visualize  crack  propagation:  a  normal 
stress  component  PuZ/xq  =  Pxx/f^o  b  shear  stress  component 
Pulido  =  Pxy/i^o  [referred  here  to  global  coordinates  with 
(X,  Y)  oriented  (normal,  parallel)  to  the  initial  crack  face] 


ture)  shown  in  Fig.  7  for  R\  =  51.  Such  trends  were  also 
found  to  be  nearly  identical  for  a  larger  inclusion  with 
Ri  =  l .51  (not  shown).  The  tendency  for  crack  deflec¬ 
tion  to  become  more  favorable  as  F /Fq  increases  is 
understandable  since  surface  energy  increases  linearly 
with  this  ratio  for  a  planar  crack  splitting  the  inclu¬ 
sion.  Regarding  the  effect  of  k/kp,  phase  held  predic¬ 
tions  in  Fig.  7  agree  with  trends  derived  analytically  via 
linear  elastic  fracture  mechanics  (He  and  Hutchinson 
1989)  for  penetration  or  deflection  of  a  straight  crack 
by  a  planar  interface  separating  two  isotropic  linear 
elastic  materials:  deflection  becomes  more  favorable 
as  the  modulus  of  the  potentially  penetrated  material 
increases. 


5  Conclusions 

A  new  phase  held  theory  of  fracture  has  been  imple¬ 
mented  in  3D  simulations,  with  results  validated  for 
isotropic  and  anisotropic  cracking  in  a  homogeneous 
material  under  distinct  mode  I  and  mode  II  far-fleld 
loading  conditions.  Potential  crack  deflection  around  a 
second-phase  inclusion  has  been  studied,  demostrating 
effects  of  relative  stiffness  and  strength  of  the  inclu¬ 
sion  to  that  of  the  surrounding  matrix  material.  These 
results  are  thought  profound  with  regards  to  the  possi¬ 
ble  design  of  tailored  nanocomposites  with  increased 
macroscopic  toughness  associated  with  crack  bridg¬ 
ing:  second-phase  particles  should  be  selected  with  the 
appropriate  properties  to  induce  deflection,  noting  that 
a  trade-off  usually  exists  between  stiffness  and  tough- 
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Fig.  6  Phase  fleld  prediction  of  crack  attraction  a  towards  soft 
(k/ko  =  0.5)  second-phase  inclusion  and  shielding  b  by  stiff 
(k/ko  =  10)  second-phase  inclusion.  In  each  case,  the  same  far- 
fleld  mode  I  boundary  conditions  are  applied,  R\  =  51,  F /  Fq  = 
1,  and  failed  material  {rj  >  0.7)  is  removed  to  visualize  crack 
propagation 


ness  of  candidate  particulate  materials  (e.g.,  stiff  but 
weak  ceramic  inclusions  versus  compliant  yet  tough 
polymer  inclusions).  In  other  words,  the  present  study 
of  crack  penetration  versus  deflection  through  inclu¬ 
sions  illustrates  how  the  phase  fleld  method  might  be 
applied  to  gain  understanding  of  how  nano-composites 
may  be  tailored  for  improved  overall  toughness.  Future 
work  shall  consider  additional  geometries  and  prop¬ 
erty  sets,  the  latter  allowing  for  interfacial  decohesion 
(Evans  1974;  Xu  and  Needleman  1993)  and  cleavage 
anisotropy  (Schultz  et  al.  1994),  both  of  which  can  be 
accommodated  by  the  present  theory. 
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Fig.  7  Predicted  tendency  for  mode  I  crack  to  penetrate  (contour 
value  0)  versus  deflect  (contour  value  1)  around  spherical  inclu¬ 
sion  (radius  R\  =  51)  with  fracture  energy  F  and  initial  stiffness 
k,  where  /q  and  Xq  are  corresponding  fixed  properties  of  the 
matrix.  Partial  defiection  is  denoted  by  contour  values  between 
0  and  1 

The  phase  field  method  is  not  restricted  to  simple 
geometries  and  elastic  linearity  typical  of  known  ana¬ 
lytical  solutions;  crack  propagation  is  insensitive  to 
mesh  construction;  and  no  spurious  parameters  control¬ 
ling  crack  propagation  direction  or  opening  displace¬ 
ment  are  required.  Since  ceramic  crystals  of  interest 
such  as  alumina  and  silicon  carbide  exhibit  competition 
among  fracture,  twinning,  and  glide  of  partial  disloca¬ 
tions  (Clayton  2009,  2010),  future  work  will  address 
such  mechanisms  simultaneously  by  assigning  multi¬ 
ple  order  parameters,  effectively  merging  the  present 
theory  with  that  in  Clayton  and  Knap  (201  la). 
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